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ABSTRACT 

We monitor the star HD 149026 and its Saturn-mass planet at 8.0 /xm over slightly more than half 
an orbit using the Infrared Array Camera (IRAC) on the Spitzer Space Telescope. We find an increase 
of 0.0227% ± 0.0066% (3.4(t significance) in the combined planet-star flux during this interval. The 
minimum fiux from the planet is 45%± 19% of the maximum planet flux, corresponding to a difference 
in brightness temperature of 480 ± 140 K between the two hemispheres. We derive a new secondary 
eclipse depth of 0.0411%±0.0076% in this band, corresponding to a dayside brightness temperature of 
1440 ± 150 K. Our new secondary eclipse depth is half that of a previous measurement (3. Oct difference) 
in this same bandpass by Harrington et al. (2007). We re-fit the Harrington et al. (2007) data and 
obtain a comparably good fit with a smaller eclipse depth that is consistent with our new value. In 
contrast to earlier claims, our new eclipse depth suggests that this planet's dayside emission spectrum 
is relatively cool, with an 8 /im brightness temperature that is less than the maximum planet-wide 
equilibrium temperature. We measure the interval between the transit and secondary eclipse and find 
that that the secondary eclipse occurs 20.9lg5 minutes earlier (2.9cr) than predicted for a circular 

orbit, a marginally significant result. This corresponds to ecos(w) — — 0.0079tQ QQ25 where e is the 
planet's orbital eccentricity and uj is the argument of pericenter. 

Subject headings: infrared: techniques: photometric - eclipses - stars:individual: HD 149026b - plane- 
tary systems 



1. INTRODUCTION 

The planet orbiting HD 149026 is unique among the 
ranks of transiting extrasolar planets. It has a mass 
comparable to that of Saturn (Sato et al. 2005; Winn 
et al. 2008; Carter et al. 2009) but its small radius and 
correspondingly high average density (Sato et al. 2005; 
Charbonneau et al. 2006; Winn et al. 2008; Nutzman et 
al. 2009; Carter et al. 2009) suggest that, unlike Saturn, 
an incredible 50 — 90% of this planet's mass must exist 
in the form of a solid icy or rocky core (Sato et al. 2005; 
Fortney et al. 2006; Ikoma et al. 2006; Broeg & Wuchterl 
2007; Burrows et al. 2007a). Together with the Neptune- 
mass planets GJ 436b (Butler et al. 2004; Gillon et al. 
2007b) and HAT-P-llb (Bakos et al. 2009), this places 
HD 149026b in a class that is distinct from that of the 
more massive "hot Jupiter" transiting planets with their 
primarily hydrogen-helium compositions. 

In light of its large rock or ice core, it is possible that 
HD 149026b also has an atmospheric composition that 
differs significantly from those of the hot Jupiters. In 
the solar system, there is a correlation between atmo- 
spheric metallicity and the percentage of planet mass 
that is core (e.g. Fodders 2003). For example, Uranus 
and Neptune have a C/H ratio of 30 — 40 times solar 
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while Jupiter's is roughly three times solar. If this re- 
lation holds true for extrasolar planets, it would sug- 
gest that HD 149026b may have a significantly metal- 
enriched atmosphere. The HD 149026 primary is also 
metal-enriched ([Fe/H]= 0.36 ± 0.05, Sato et al. 2005), 
thus even if HD 149026b's atmosphere simply refiects the 
composition of its primordial nebula we would expect it 
to have a slightly higher heavy metal content than a typ- 
ical hot Jupiter. 

By observing the decrease in light as the planet passes 
behind its parent star in an event known as a sec- 
ondary eclipse, we can characterize the properties of 
HD 149026b's dayside emission spectrum. Harrington 
et al. (2007) (hereafter H07) measured HD 149026b's 
secondary eclipse at 8 /im and found a depth of 
0.084%^[J p°2%j corresponding to a brightness tempera- 
ture^ of 2300 ± 200 K for the planet. This is significantly 
higher than expected, as HD 149026b has an equilibrium 
temperature^ of only 1740 K (H07). In order to repro- 
duce this fiux in the case where the planet emits as a 
blackbody, we must assume that the planet absorbs and 
then instantaneously re-radiates all of the incident flux 
from its star on the day side. In this scenario the planet's 
night side would have effectively zero flux; this follows 
from the requirement that the energy emitted by the 
planet does not exceed the energy it absorbs (note that 
residual heat from formation is not expected to play a 

^ Brightness temperature is defined in this case as the temper- 
ature required to match the observed planet-star flux ratio in the 
8 fim Spitzer IRAC bandpass assuming that the planet radiates as 
a blackbody and using a Kurucz atmosphere model for the star. 

* Equilibrium temperature is calculated by assuming that the 
planet absorbs all incident flux (i.e., zero albedo) and then re- 
radiates that energy over its entire surface as a blackbody of the 
stated temperature. 
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significant role in the energy budgets of elose-in planets). 
Tliis is not unreasonal^le. as tlie time scale for tidal lock- 
ing of close-in planets such as HD 149026b is consider- 
ably shorter than the ages of these systems (Guillot et al. 
1996; Lubow et al. 1997). Uneven irradiation of the up- 
per atmosphere could produce a large observed day-night 
temperature gradient if radiative time scales are much 
less than advective time scales at the level of the mid-IR 
photosphere (i.e., flux is absorbed and re-emitted by gas 
on the dayside in less time than it takes that same parcel 
of gas to reach the nightside hemisphere). If the planet's 
spectrum differs significantly from that of a blackbody 
with enhanced emission at 8 /im, the constraints on the 
night-side emission becomes correspondingly less strin- 
gent. If we compare HD 149026b's 8 /xm flux to the 
predictions from ID atmosphere models for the planet 
(e.g. Fortney et al. 2006; Burrows et al. 2GG8a), we find 
that models with a temperature inversion and water fea- 
tures in emission at 8 /im still require a hot day side and 
correspondingly large day-night temperature gradient in 
order to match the 8 /zm flux observed by H07. 

In this paper we characterize HD 149G26b's 8 /im phase 
variation and corresponding day-night temperature gra- 
dient by monitoring the system continuously over slightly 
more than half an orbit, beginning before the transit and 
ending after the secondary eclipse. By measuring the in- 
crease in brightness as the planet's dayside rotates into 
view, we can estimate the day-night temperature gradi- 
ent and constrain corresponding atmospheric circulation 
models for the planet. We have previously published sim- 
ilar observations of HD 189733b at 8 and 24 /im (Knut- 
son et al. 2007, 2009a) where we invert the observed 
phase variation to produce a longitudinal temperature 
profile for the planet. These observations revealed that 
HD 189733b has a warm night side (approximately 200 K 
cooler than the 1250 K day side), but lower-cadence ob- 
servations of the non-transiting planets v And b (Har- 
rington et al. 2006) and HD 179949b (Cowan et al. 2007) 
suggest that other hot Jupiters may have large day-night 
temperature gradients. HD 149026b's lower mass and in- 
creased core fraction make it qualitatively different than 
any of these systems, and there are no published gen- 
eral circulation models for this planet analogous to those 
available for HD 189733b and HD 209458b (e.g. Show- 
man et al. 2009). As a result, we have no a priori pre- 
dictions for the nature of the day-night circulation on 
HD 149026b other than the constraint provided by the 
secondary eclipse measurement of H07. 

These same data also allow us to sc;ai'ch for time- 
varying properties of the system by comparing the depths 
and relative times of our transit and secondary eclipse to 
previously published 8 /xm Spitzer observations of the 
planet's transit (Nutzman et al. 2009, obtained in 2007) 
and secondary eclipse (H07, obtained in 2005). If the 
planet's orbit is changing over time, perhaps as the result 
of interactions with an unknown second planet in the sys- 
tem, it could cause deviations in the timing of successive 
transits (Miralda-Escude 2002; Holman & Murray 2005; 
Agol et al. 2005). Comparing our 8 /im planet-star ra- 
dius ratio to the value obtained by Nutzman et al. (2009) 
allows us to search for changes in the effective radius of 
the planet at 8 fim and, if the two values are consistent, 
obtain an improved estimate for this quantity. By mea- 
suring the interval between the transit and secondary 



eclipse we can determine a value for eco,s((x-') where e is 
the planet's orbital eccentricity and uj is the argument 
of pericenter, and compare it to the corresponding value 
from H07. Finally, by comparing the 8 /im secondary 
eclipse depth from our observations to that of H07, we 
can search for variations in the planet's 8 /xm dayside 
flux. 

In §2 we describe our treatment of the data, including 
our initial photometry, the use of a new preflash tech- 
nique to remove the detector ramp, and our fits to the 
transit, secondary eclipse, and phase curve data. In §3 
we discuss the results of these fits and check for variabil- 
ity in the relative depth and timing of our transit and 
secondary eclipse as compared to previously published 
data. We also compare our day- and night-side fluxes to 
the predictions from ID atmosphere models and calcu- 
late an energy budget for the planet. In §4 we summarize 
our results and discuss the potential impact of upcoming 
observations of this planet's secondary eclipse and phase 
curve at shorter wavelengths. 

2. OBSERVATIONS 

We monitored HD 149026b continuously for 41 hours 
using 8 /xm IRAC subarray (Fazio et al. 2004) on the 
Spitzer Space Telescope (Werner et al. 2004), from UT 
2008 May 11 — 12. Our observations began 2.5 hours 
before the start of the transit and finished 44 minutes 
after the end of the secondary eclipse. We used an inte- 
gration time of 0.4 s, well below saturation for this star, 
and obtained a total of 344,704 images. 

2.1. Initial Photometry 

Because HD 149026 is a bright star {K = 6.819, one- 
third the brightness of HD 189733 in the infrared) rela- 
tive to the background at these wavelengths, we calculate 
the fiux from the star in each image using aperture pho- 
tometry with a radius of 3.5 pixels. We determine the po- 
sition of the star in each image as the position- weighted 
sum of the flux in a circular region with a radius of five 
pixels (rounded to the nearest pixel) centered on the ap- 
proximate position of the star. We calculate position es- 
timates with larger and smaller apertures, but find that 
a five pixel aperture produced the lowest RMS in the re- 
sulting time series. Similarly, we find that decreasing our 
photometric apertures to values smaller than 3.5 pixels 
led to position-dependent fiux losses, while the time se- 
ries from larger apertures become increasingly noisy. We 
repeat our analysis for apertures ranging from 3.5 — 7 
pixels and find that our results for the measured sec- 
ondary eclipse depth and phase variation amplitude are 
consistent in all cases. We also try fitting a point spread 
function (psf) to the images, using either the in-fiight 
point response functions generated from calibration test 
data ^ or a psf derived from our own data, and conclude 
that aperture photometry produces the optimal results. 

Because the subarray is only 32 x 32 pixels, care is 
needed in the choice of regions for background flux esti- 
mates. We exclude all pixels within a circular region with 
a radius of 9 pixels centered on the position of the star, 
as well as the top row of the array which has consistently 
low flux values, and trim 3a outliers from the remaining 
set of pixels. We then make a histogram of the fluxes in 

^ Available at http://ssc.spitzer.caltech.edu/irac/psf.html 
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Fig. 1. — This figure shows the RMS of the time series be- 
tween the end of the transit and the start of the secondary eclipse 
(10 < dt < 35 hours) over a range of bin sizes (black line, bins up 
to 50,000 points). We have fitted and removed a linear function 
of time from these data to correct for the effects of the planet's 
phase variation before binning the data and calculating the result- 
ing RMS. The solid red line shows the RMS of the unlainned data 
scaled by y'n where n is the number of points in each bin. The 
increased noise for bin sizes between 1000 — 4000 points (10 — 30 
minute time scales) may be related to the pointing jitter of the 
telescope, which shifts back and forth by approximately 0.1 pixels 
in X and y over 1 hour cycles. Although we find no evidence for a 
direct relation between the measured x/y positions of the star on 
the array and the corresponding stellar fluxes, the fact that this 
excess noise disappears for bin sizes larger than 8,000 points (1 
hour intervals) suggests these two phenomena may be related. 

the remaining pixels and fit a Gaussian function to the 
central 7 of 10 bins in this distribution (further reduc- 
ing the effects of any remaining high- or low-flux outliers 
in the distribution) to estimate the mean background 
level in each individual image. We obtain identical re- 
sults if we include the top row of pixels, increase the 
size of the exclusion region in the center of the array, or 
exclude an additional four rows from y = 11 — 16 corre- 
sponding to particularly bright diffraction spikes on the 
star's psf. We find that the background flux during our 
observations is consistently below zero, with a median 
value of —0.053 MJy Sr~^, suggesting that the sky dark 
used for these observations is over-subtracting the back- 
ground, but this has no effect on our final photometry 
provided that the negative background value is correctly 
subtracted from each image before estimating the flux 
from HD 149026. 

We calculate the JD value for each image as the time 
at mid-exposure, and apply a correction to convert these 
JD values to the appropriate BJD, taking into account 
Spitzer's orbital position at each point during the obser- 
vations. The background and stellar fluxes in the first 10 
images and the 58th image in each set of 64 have consis- 
tently lower values than the rest of the exposures, but we 
find that this behavior disappears after the background 
is subtracted. As a check we exclude these low frames 
and repeat our analysis, and find that we obtain consis- 
tent results (secondary eclipse depth and phase variation 
amplitude) in both cases. We include these frames in our 
final analysis. 

Our final time series has a number of high-flux outliers 
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Fig. 2. — Time series for a 3.5 pixel region centered on the 
approximate position of HD 149026 (corresponding to the center 
of the subarray) over the 30 minutes of our preflash observations of 
M17. The average flux in this region is 7400 MJy Sr~^, but we see 
an identical behavior for regions near the corners of the array with 
average fluxes ranging from 3900 — 8200 MJy Sr~^. This suggests 
that we are effectively saturating out the detector ramp in the first 
five minutes of our preflash observations. 

produced by the presence of transient hot pixels in the ar- 
ray. We remove these points by discarding images where 
the X position, y position, or total flux from the star in 
a 3.5-pixel aperture differed by more than 3cr from the 
median value for these quantities in each set of 64 images 
(because we calculate the star's position as the position- 
weighted sum of the flux in a given region, transient hot 
pixels may affect these estimates even if they are not in- 
cluded in the final photometry). We trim 0.78% of the 
images using this method, and see no evidence for re- 
maining outliers in the final time series. We calculate 
the RMS of the points in each bins in the phase curve 
plotted in Fig. 7 and find that the average RMS in these 
bins is 0.826%, which is a factor of 1.27 higher than the 
predicted photon noise from the star alone (we do not 
include the noise contribution from the background flux 
because this flux is negative in our images, and we do 
not have a reliable estimate of its true value with a cor- 
rect sky dark subtraction). This is consistent with pre- 
vious Spitzer photometry in this bandpass for GJ 436, 
HD 149026, and HD 189733 (Deming et al. 2007; Har- 
rington et al. 2007; Knutson et al. 2007), all of which have 
roughly comparable fluxes in this bandpass. We find that 
this noise is reasonably white; Fig. 1 shows the RMS of 
the time series between the end of the transit and the 
start of the secondary eclipse (10 < c?t < 35 hours) over 
a range of bin sizes up to 50,000 points after the data 
has fitted with a linear function of time which was then 
removed to ensure that the real astrophysical change in 
the planet's flux does not inflate the resulting RMS esti- 
mates. 

2.2. Preflashing and the Detector Ramp Effect 

After determining the optimal method for estimat- 
ing the fluxes in each bandpass, we must remove any 
remaining detector effects from the data. There is a 
well-documented detector effect (Harrington et al. 2007; 
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Knutson et al. 2007, 2008; Charbonneau et al. 2008) 
in this array that causes the effective gain (and thus 
the measured flux) in individual pixels to increase over 
time. This effect has been referred to as the "detector 
ramp", and has also been observed to occur in the IRS 
16 fim peak-up and MIPS 24 fim arrays, which are made 
from the same material (Deming et al. 2006; Charbon- 
neau et al. 2008). The IRAC 5.8 /im array also shows 
a related behavior, with an initial upward asymptote in 
the measured flux over the first 30 — 60 minutes of ob- 
servation followed by a more gradual downward slope 
(Charbonneau et al. 2008; Knutson et al. 2008, 2009b; 
Machalek et al. 2008). The size of this effect depends on 
the illumination level of the individual pixel. Pixels with 
high illuminations (>250 MJy Sr^^ in the 8 /im channel) 
will converge to a constant value within the first hour of 
observations, whereas lower-illumination pixels will show 
a linear increase in the measured flux over time with a 
slope that varies inversely with the logarithm of the il- 
lumination level. It has been suggested that this effect 
may be the result of charge-trapping in the array, with 
higher-illumination pixels filling their wells more quickly 
than lower-illumination pixels. 

We mitigated this effect by observing M17, a nearby 
open cluster with diffuse bright emission at 8 /im, 
centered on the position a = 18''20™28.0^ S = 
— 16°12'19.5" for 32 minutes immediately prior to the 
start of our science observations. This preflash is de- 
signed to saturate out the ramp effect by exposing the 
array to a uniformly bright source with fluxes that are 
approximately 10 x higher than the peak flux in the cen- 
ter of HD 149026's psf. The region of M17 that we use for 
our preflash has fluxes between 3500 — 8000 MJy Sr~^ 
in the IRAC 8 /xm subarray aperture, with an average 
flux of 7400 MJy Sr~^ in a circular region with a radius 
of 3.5 pixels centered on the position of HD 149026 in 
our science images. We use the same exposure time and 
subarray observing mode for our preflash that we use for 
our observations of HD 149026, obtaining a total of 4480 
preflash images. The peak flux in these images occurs at 
the position of a star in the corner of the array with a 
maximum flux of 13,000 MJy Sr^^ or 100,000 e", cor- 
responding to half of the full well depth in this array. 
Because we are well below saturation, we do not expect 
any latency effects in our subsequent science images. 

We evaluate the effectiveness of our preflash by calcu- 
lating the total flux in an aperture with a 3.5 pixel radius 
centered on (1) the center of the array in the preflash im- 
ages, corresponding to the position of HD 149026 in our 
science observations, (2) a region near the corner of the 
array in the preflash images with a lower flux around 
3500 MJy Sr~^, and (3) a bright star in the opposite 
corner of the array in the preflash images with a peak 
flux of 13,000 MJy Sr'^. We find that in all cases there 
is a linear increase in the measured flux during the first 
five minutes of observations with a total amplitude of 
0.4% (see Figure 2), followed by a slight decrease in the 
measured flux and then remaining effectively constant for 
the final 15 minutes of observations. Based on these data 
we conclude that the effectiveness of our preflash is inde- 
pendent of the brightness of our preflash source above a 
minimum threshold value of 4000 MJy Sr~\ and that the 
charge traps across our entire 32 x 32 pixel subarray are 
flUed in the first five minutes of our preflash observations. 
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Fig. 3. — Final time series for HD 149026 (black points) over the 
41 hours of our observations, with data binned in 428 s intervals 
containing approximately 1000 points per bin. There is still a small 
ramp at early times, but its amplitude is smaller than the equiva- 
lent ramp with no preflash and the ramp converges to a constant 
value within the first 10 hours of observations. For comparison we 
show equivalently binned photometry for a transit observation of 
HD 149026b observed on UT 2007 August 14 in this same bandpass 
with identical 0.4 s exposures in the 8 ^m IRAC array (Nutzman 
et al. 2009). Unlike our data these older data were not preflashed, 
and the increased slope of the resulting photometry indicates the 
amplitude of the unmitigated detector ramp for this star. 

When we plot our final time series for HD 149026 (Fig. 
3), we find that the preflash has moderated the ramp ef- 
fect in our science data but it has not removed it entirely. 
A comparison to identical observations of HD 149026 
(8 fim IRAC array, 0.4 s integration times, Nutzman et 
al. 2009) from 2007 reveals that the amplitude of the 
ramp in our preflashed data is considerably smaller than 
the equivalent ramp in these older data. We speculate 
that the smaller ramp in our data may be the result of a 
slow leakage of charge from the filled charge traps during 
the 14.4 minute interval between the end of our preflash 
and the start of our science data that partially resets the 
detector ramp effect. We note that unlike the original 
unmitigated ramp, we see no variation in the slope of 
this new ramp with increasing aperture, and the back- 
ground fluxes show a slight decrease in flux during the 
flrst few hours of observations. 

Our data do not show any evidence for a detector ramp 
after the first 10 hours of observations (see §2.4 for a 
more detailed discussion on this point). We account for 
the ramp at early times (< 10 hours after the start of 
observations) by fitting the following quadratic function 
of In (dt) simultaneously with our transit fit: 

/ = ci + C2 In {dt + a) + C3 (In {dt + Ci)f (1) 

where / is the measured flux, dt is the elapsed time in 
days since the start of the observations, and ci — C4 are 
free parameters in the fit. We also tried fitting this slope 
with a linear function of In {dt) or a simple exponential 
function, but found that Eq. 1 provided the best fit to 
the data. 

2.3. Eclipse Fits 
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Fig. 4. — This plot shows the first ten hours of data with the 
detector ramp removed and the best-fit transit light curve over- 
plotted In this plot we use the fit with o/il« and i fixed to their 
best-fit values from Carter et al. (2009). The data are binned in 
259 s intervals containing approximately 600 points per bin. The 
uncertainties on each individual bin arc calculated as the standard 
deviation of the flux values in that bin divided by the square root of 
the number of points. Data during the center of transit show sev- 
eral above-average deviations from the best-fit transit light curve; 
the time scale of these deviations is similar to the planet crossing 
time scale, and such deviations could be caused by the presence of 
star spots intersecting the planet's path across the surface of the 
star. 

We fit both Eq. 1 and the transit curve to the first 

ten hours of data simultaneously using a Markov Chain 
Monte Carlo (MCMC) method (e.g. Tegmark et al. 2004; 
Ford 2005; Winn et al. 2007a; Bnrkc et al. 2007). We ex- 
clude the first 30 minutes of data from our fit, as the ramp 
is particularly steep during this time, and set the uncer- 
tainties on individual points equal to the RMS variance 
of the time series between 0.8 — 1.0 days after the start 
of the observations (this is the middle of the time series, 
well after the transit but before the start of the secondary 
eclipse, and we have selected a shorter segment to avoid 
inflated errors due to the planet's phase variation). We 
fix the four parameter nonlinear limb-darkening coeffi- 
cients (Clarcit 2000) to the values used by Nutzman et 
al. (2009), which are derived by fitting the appropriate 
limb-darkening law to a Kurucz stellar atmosphere model 
for a Teff = 6250 K, log{g) = 4.5,[Fe/H]= 0.3 star (Ku- 
rucz 1979, 1994). We calculate our transit fight curve 
using the routines from Mandel & Agol (2002) and fit for 
four free parameters in addition to the four parameters 
in Eq. 1, including: the planet-star radius ratio Rp/R^, 
the ratio of the planet's orbital distance to the stellar ra- 
dius a/R^, the orbital inclination z, and the transit time. 
These four parameters are independent of assumptions 
about the properties of the star, and can be derived di- 
rectly from the light curve. Using Kepler's law together 
with the planet's orbital period and the stellar mass, we 
can then derive the corresponding values for the stellar 
and planetary radii ii, and Rp. This gives us a total of 



Fig. 5. — This plot shows the 7.2 hours of data used for the 
secondary eclipse fit, beginning 0.2 days from the predicted center 
of the eclipse). The data are binned in 517 s intervals containing 
approximately 1200 points per bin. The uncertainties on e3x:h. in- 
dividual bin are calculated as the standard deviation of the flux 
values in that bin divided by the square root of the number of 
points. The two overplotted curves show the best-fit eclipse light 
curve with best-fit transit time (solid line) or with the time fixed 
to its predicted value (dashed line). We find that the secondary 
eclipse occurs 20.9^g g minutes earlier than predicted for a circular 
orbit (see discussion in §3.1.2). 

eight free parameters in our fit to the transit, and we run 
our Markov Chain for a total of 10^ steps. 

After running the chain, we search for the point in the 
chain where the value first falls below the median of 
all the values in the chain (i.e. where the code had 
first found the best-fit solution), and discard all the steps 
up to that point. We take the median of the remaining 
distribution as our best-fit parameter, with errors calcu- 
lated as the symmetric range about the median contain- 
ing 68% of the points in the distribution. We find that 
a/i?* and i are highly correlated and have asymmetric 
histograms, therefore for these two parameters we calcu- 
late separate positive and negative error bars. 

As noted above, our transit data are virtually identical 
to that of Nutzman et al. (2009) with the sole exception 
of a fiatter detector ramp (see Fig. 3). Carter et al. 
(2009) recently presented a second, independent analysis 
of the Nutzman et al. (2009) data in combination with 
new observations of the transit between 1.1 — 2.0 fiia us- 
ing NICMOS on the Hubble Space Telescope. The au- 
thors in this study derive improved constraints on a/i?* 
and i relative to an analysis of the Spitzer data alone. In 
order to facilitate comparisons between these two data 
sets, we repeat our fits with these two parameters fixed 
to their best-fit values from Carter et al. (2009) and allow 
the remaining six parameters (four from Eq. 1, Rp/R^, 
and the transit time) to vary freely. We give the results 
of both fits in Table 1 below. 

We perform our fit to the secondary eclipse using sim- 
ilar methods; there is no detector ramp in this region 
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TABLE 1 

Best-Fit Transit and Secondary Eclipse Parameters 



Parameter 



Transit Parameters 

-Rp/i?* 

i n 

Center of Transit 
O-C (minutes)'' 

Secondary Eclipse Parameters 
Depth 

Center of Eclipse 
O-C (minutes)'i'« 



Independent Fit^ 



0.05216 ± 0.00078 
88.0lJ-| 

2454597.70721 ± 0.00040 BJD 
1.10 ±0.63 



0.0411% ± 0.0076% 
2454599.1304 ± 0.0048 BJD 

-20.91H 



Carter et al. (2009) Fit'' 



0.05253 ± 0.00076 
84.50 
5.99 

2454597.70709 ± 0.00039 BJD 
0.92 ± 0.62 



^Thcsc values are from an independent fit to our data alone, with no prior assumptions about a/R^ or i. 
^Hcrc we set fix a/R^ and i to their best-fit values from Carter ct al. (2009). 

^The probability distribution for a/R-j, is highly asymmetrie, so we use the mode instead of the median 
as our best-fit value and set our uncertainties to the 16th and 84th percentiles of the distribution. 

^'Calculated using the new ephemeris derived in this paper; the uncertainties include the uncertainties in 
both the observed and calculated transit times. 

^The O-C for the secondary eclipse is calculated assuming that the secondary eclipse occurs half an orbit 
after the transit plus an additional 45 s delay from the increased light travel time in the system (Loeb 
2005; Harrington et al. 2007). 



of the time series so we limit our out-of-transit normal- 
ization to a simple linear function of time which we fit 
simultaneously with the eclipse. We select the final 7.2 
hours of data used for the secondary eclipse fit (begin- 
ning 0.2 days before the predicted center of the eclipse 
in order to provide an adequate out-of-eclipse baseline). 
We find the slope of the linear term in this fit is con- 
sistent with zero, and obtain similar results if we fit a 
constant value to the out-of-eclipse data (see discussion 
in Kmitson ct al. 2009a). We allow both the depth and 
timing of the secondary eclipse to vary independently, 
and take the other parameters for the system (planetary 
and stellar radii, orbital period, etc.) from Nutzman et 
al. (2009). Given the relatively low signal-to-noise of this 
eclipse as compared to the transit, we conclude that up- 
dating the system parameters to the values from Carter 
et al. (2009) would have no eff'ect on our final best-fit 
transit depth and corresponding uncertainties. As be- 
fore, we calculate our eclipse curve using the equations 
from Mandel & Agol (2002) with limb-darkening set to 
zero. 

Figures 4 and 5 show the final binned data from our 
fits to the transit and secondary eclipse with the out- 
of-transit trends removed and the best-fit eclipse curves 
ovcrplottcd. The corresponding best-fit parameters are 

given in Table 1. 

2.3.1. Fits to 8 iim Secondary Eclipse Data from 2005 

As part of this paper we would like to compare our 
new 8 /im secondary eclipse depth to the 8 fim secondary 
eclipse depth from 2005 presented in H07. Over the past 
two years there have been a number of changes both to 
the standard Spitzer pipeline and also in the accepted 
best-practice methods for dealing with these data. In 
light of these changes, we decided to evaluate the robust- 
ness of the measured 2005 eclipse depth by downloading 
the images^° from 2005 and performing our own indepen- 
dent analysis of these data. The images we downloaded 

'" These data are publicly available through the Spitzer archives. 



from the Spitzer archives, which are periodically repro- 
cessed with updated versions of the Spitzer pipeline, were 
last updated in June 2006 (version S14.0.0). We use these 
images to create an updated photometric time series, fol- 
lowing the same methods as those described in §2.1, and 
compare the best-fit eclipse depths obtained using these 
data with those obtained using the photometry from the 
original H07 paper'^'^. We find that we obtain the same 
eclipse depths using both our new photometry and the 
original H07 photometry. As described above, we calcu- 
late our aperture photometry using a circle with a radius 
of 3.5 pixels, but we obtain consistent results using the 
slightly larger 4.0 pixel aperture selected by H07. 

Unlike our data, the H07 data cycles between nine dif- 
ferent nod positions, with a total of twelve cycles through 
these nod positions (see the supplemental methods sec- 
tion of H07 for more information). This creates a rather 
complicated time series, as there is both an overall ramp 
in the data comparable to the one described in §2.2, and 
a series of smaller ramps that begin each time the tele- 
scope moves to a new pointing position. We follow the 
general method described in H07 and fit these effects 
with a combination of functions, including a function to 
describe the overall ramp: 

/™„^p = ci(l-e-('^*-+°-^5)) (2) 

where ci and C2 are free parameters in the fit, dtobs 
is the elapsed time from the start of the observations in 
days, and 0.75 is an arbitrary time offset. We tried fits 
in which this time offset was a free parameter, but the 
additional degree of freedom did not improve the qual- 
ity of the resulting fit and our Markov chains indicated 
that the time offset and C2 were highly degenerate. We 
selected the value of 0.75 from the center of the distribu- 
tion of values where this parameter was allowed to vary 
freely. We also ran fits using a quadratic functions of 

" This photometry was provided as part of the H07 Nature 
paper's supplemental information. 
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In dt, similar to that given in Eq. 1 but with the phase 
shift C2 in that equation fixed to a value of 0.02 days, 
and found that the exponential function was a better fit 
to the data (H07 reached a similar conclusion). 

As mentioned above, there is a series of smaller ramps 
in these data that reset each time the telescope moves 
to a new pointing position. We fit these ramps with a 
function of In. dt, comparable to that of Eq. 1: 



nod 



I + C3 In {dtnod + C4) 



(3) 



where C3 and C4 are free parameters and dtnod is the 
elapsed time in days since the start of observations at the 
current nod position (this parameter resets to zero each 
time the telescope repoints to the next nod position in 
the cycle). Our choice to fit this ramp with an analj^tic 
function differs from that of H07, who defined the ramp 
at each nod position using spline interpolation through 
the photometry from twelve images in each visit where 
the final three observations were fixed to unity (see §5 in 
the supplemental information section in H07). 

We also include constant flux offsets C5 — C12 at eight of 
the nine nod positions (the offset at the first nod position 
is fixed to zero) as free parameters in our fits in order to 
account for position-dependent differences in the mea- 
sured fiuxes across the array. H07 do the same in their 
fits. We allow both the depth and timing of the sec- 
ondary eclipse to vary independently, and take the other 
parameters for the system (planetary and stellar radii, 
orbital period, etc.) from Nutzman et al. (2009) as de- 
scribed in §2.3. This choice differs from that of H07, who 
describe their eclipse with four free parameters (eclipse 
duration between half-light points, center time, depth, 
and limb-crossing duration), as the uncertainties in the 
eclipse duration and length of ingress/egress were much 
higher at the time this paper was written. 

We carry out our fits using a Markov Chain Monte 
Carlo method, identical to the one described in §2.3 with 
a total of 14 free parameters (H07 have 23 free parame- 
ters in their fit): 

f measured ~ frampipl^ ^2) * fnodi^^} C4) (4) 
*/ecHpse(ci3,Cl4) + foffset{c5 - C12) 

where fmeasured IS the measured flux, framp is given 
by Eq. 2 and describes the overall ramp in the data, 
fnod is given by Eq. 3 and describes the smaller ramp at 
each new nod position, /offset is the constant flux off- 
set for the last 8 of 9 total nod positions, and f eclipse is 
the secondary eclipse function (calculated the same way 
as the secondary eclipse in §2.3) where the eclipse depth 
ci3 and the time of center of eclipse C14 are left as free 
parameters. We set the uncertainty on each individual 
point equal to the same value used for our new photom- 
etry in §2.1 (25.790 MJy Sr~^ for a 3.5 pixel aperture or 
26.804 MJy Sr~^ for a 4.0 pixel aperture), as this reflects 
the actual RMS achieved for this star in the absence of 
any significant detector effects. Our MCMC fit is differ- 
ent than the fitting method employed by H07, who used 
the Levenberg-Marquardt method to minimize their 
and determine the corresponding best-fit solution, with a 
bootstrap method with 17,030 trials (as opposed to our 
10^ step Markov chain) to estimate the corresponding 
uncertainties. 
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Fig. 6. — This plot shows the updated photometry for an 
8 lan cchpsc of HD 149026b observed on UT 2005 August 24, 
originally published by Harrington et al. (2007). The time series 
shown in this plot is calculated our own photometry pipeline with 
a 4.0 pixel radius aperture and has detector effects removed. These 
observations lasted 6.0 hours, but we have plotted the same time 
interval (7.2 hours) as in Fig. 5 to facilitate comparisons. The data 
are binned in 527 s intervals containing approximately 1200 points 
per bin. The uncertainties on each individual bin arc calculated in 
the same way as Fig. 5. The two ovcrplottcd curves show the best- 
fit eclipse light curve from our analysis (solid line) and an eclipse 
curve calculated using the best-fit eclipse depth and time from H07 
(dashed line). Wc conclude that these data arc consistent with the 
secondary eclipse depth from our 2008 data, and find no evidence 
for time variability in the planet's 8 /im flux. 

As discussed in §2.1, images taken in subarray mode 
arc stored in sets of 64 and the background and stel- 
lar fluxes in the first 10 images and the 58th image in 
each set of 64 have consistently lower values than the 
rest of the exposures. For our new data we found that 
this effect was removed by the background subtraction, 
but as a check we repeat our fits to the updated pho- 
tometry for the 2005 data with and without these im- 
ages. We find that we obtain consistent eclipse depths 
when we trim the 58th image alone (an eclipse depth 
of 0.055% ± 0.007%), trim the 1st five images and 58th 
image (an eclipse depth of 0.057% ± 0.007%), or the 
1st 10 images and the 58th image (an eclipse depth of 
0.060% ± 0.008%). If we change our method for remov- 
ing transient hot pixels from the time series from the one 
described at the end of §2.1 to a simple flux cutoff (any 
images with fluxes lower than 3000 MJy Sr"-*^ or higher 
than 3230 MJy Sr^^ for a 3.5 pixel radius aperture are 
removed. whc;re the flux c;utoffs are chosen using a simple 
visual inspection of the time series), we obtain an even 

deeper eclipse depth of 0.086%+g []2o%- Using a 4.0 pixel 
aperture^ and our original bad pixel filtering method, and 
either trimming both the 58th image and 1st 10 images 
in each set of 64 or leaving those images in, we obtain 
eclipse depths of 0.052% ± 0.009% and 0.048% ± 0.008%, 
respectively. If we allow the time offset of 0.75 in Eq. 2 
to vary freely in our fits with a 3.5 pixel aperture and 



8 



1.0008 



1.0006 

X 

il 1.0004 

CD 
> 

I 1.0002 

DC 

1.0000 
0.9998 

-0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 

Orbital Phase 

Fig. 7. — This plot shows the data used in the phase curve fit, beginning 10 hours after the start of the observations (to avoid the initial 
detector ramp) and excluding data within 8 minutes (2cr uncertainty in the original best-fit eclipse time, see §3.1.2) of the start and end of 
the best-fit eclipse solution. The stellar flux as measured at the center of the secondary eclipse has been normalized to unity (dashed line). 
The data before the start of the secondary eclipse are binned in 2.2 hour intervals containing approximately 20,000 points per bin, and 
the error bars for each bin are calculated as the standard deviation of the flux values in that bin divided by the square root of the number 
of points. This is a slight over-estimate of the uncertainties on individual points, as the decreasing RMS for bin sizes greater than 1 hour 
(Fig. 1) indicates that there is a component of correlated noise, perhaps connected to the pointing jitter of the telescope, that averages out 
on longer time scales. Data after the end of the secondary eclipse are binned into a single point, and we show the best-fit secondary eclipse 
depth with la uncertainties at the center of the best-fit eclipse solution for comparison. The solid red line is the best-fit two-hemisphere 
model where the locations of the two hemispheres are allowed to vary freely in the fit, and the dotted red line is the same model with the 
locations of the two hemispheres fixed on the substellar and antistellar points. The additional degree of freedom in the first (solid line) 
model reduces the of the fit by 6.2, a clear improvement. The two blue lines are for model fits in which the day-night fiux distribution 
varied as the cosine of the longitude from the substellar point with (solid line) and without (dotted line, located underneath the red dotted 
curve) an additional phase shift as a free parameter in the fit. The two- hemisphere and sinusoidal models are indistinguishable at the 
precision of the data. 




trim the 1st 10 images and the 58th images we obtain 
a shghtly deeper eclipse depth with larger uncertainties 
(0.076%+° °20%)- However, this fit is not stable, as we 
obtain a depth of 0.058%+° °°^| if we include the 1st 
10 frames in this fit, and if we use photometry with a 
slightly larger 4.0 pixel aperture instead of 3.5 pixels we 
obtain a value of 0.052% ± 0.008% for the eclipse depth. 
We plot one example of our binned photometry (best-fit 
eclipse depth of 0.048% ± 0.008% with detector effects 
removed and best-fit eclipse curve overplotted in Fig. 6. 
The fit plotted in this figure uses photometry with a 4.0 
pixel aperture, our best bad pixel filtering method, in- 
cludes the first 10 and 58th images in our fit, and fixes 
the time offset in Eq. 2 to 0.75. 

The diversity of echpse depths (0.05% - 0.09%) ob- 
tained in these fits suggest that the final result is sensi- 
tive to our specific choice of functions, fitting routines, 
and bad pixel trimming methods. This is not surprising 
in light of the many degrees of freedom required to fit 
the detector effects in these data. As noted above, we 
obtained consistent results using both our updated ver- 
sion of the 2005 photometry (with our own photometry 
routines) and the original photometry provided by HOT. 

As a final test we compare the value for the same 
kind of fit plotted in Fig. 6, but using the original HO 7 
photometry instead of our new photometry, to the 
value of the best-fit solution obtained by H07. We trim 
the same images from our final time series as H07; the 
only difference is that we set the uncertainty on individ- 



ual points equal to the average uncertainty from H07, 
rather than using the vector of individual uncertainties 
provided by H07. Our best-fit solution as determined 
from a 10^ step Markov chain, has an eclipse depth of 
0.050% ± 0.008% and a reduced value of 0.67359847 
with 48, 145 points in the fit and 14 degrees of freedom. 
This is slightly smaller than the reduced x^ of 0.67368895 
obtained by H07 with the same number of points and 23 
degrees of freedom (J. Harrington 2009, private commu- 
nication), but the difference is not statistically signifi- 
cant. 

Based on this analysis, we conclude that these data 
are consistent with a shallower eclipse depth than that 
published by H07. As we discuss in more detail in §3.1.1, 
we find that there is no evidence for variability in the 
planet's 8 /im flux over the three- year interval between 
these data and our new observations. 

2.4. Phase Curve Fits 

Due to the small amplitude of the secondary eclipse 
depth (0.041% ±0.008%), a robust detection of the phase 
variation of HD 149206b in these data presents a con- 
siderable challenge. For comparison, the measured am- 
plitude of the phase variation for HD 189733b, whose 
host star is three times brighter than HD 149026 in this 
bandpass, was 0.12% ± 0.02% (Knutson et al. 2007). 
This is five times larger than the best-fit amplitude for 
HD 149026b's phase variation given below. 

As discussed in §2.2, we trim the first 10 hours of data 
from our phase curve fits, removing the part of the time 
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TABLE 2 

Minimum and Maximum Hemisphere- Averaged Fluxes 



Parameter 



Planet-Star Flux Ratio Brightness Temperature'' 



Linear Fit to Phase Curve 

(X^ = 228872.7, 228920 points & 2 DOF) 

F b 



0.0411% ± 0.0076% 
0.0185% ± 0.0100% 
45% ± 19% 



1440 ± 150 K 
960 ± 240 K 

480 ± 140 K 



Two Hemisphere Model Fit 

(X^ = 228867.6, 228920 points & 3 DOF) 



0.0539% ± 0.0138% 
0.0161% ± 0.0080% 
30% ± 13% 



1680 ± 260 K 
910 ± 200 K 

780 ± 240 K 



^Here we assume the planet emits as a blackbody and solve for the temperature that reproduces the 
observed planet-star flux ratio in our 8 /Ltm bandpass. We use a PHOENIX NextGen model atmosphere 
(e.g. Hauschildt et al. 1999) for the star with an effective temperature of 6160 K and log(g)= 4.19; this is 
the same model used to calculate the planet-star flux ratios for the ID models in Figures 11 and 12. 

^The maximum hemisphere-averaged flux (given as a percentage of the stellar flux) in this fit is set equal 
to the dayside flux as determine from the secondary eclipse depth. 

*^The uncertainty in the minimum hemisphere-averaged flux includes the uncertainties from both the 
secondary eclipse depth estimate and the relative increase in flux (in this case, the slope of the line in the 
phase curve flt). 

^Thc maximum flux in this flt occurs after the end of our observations, and the uncertainties in its best-flt 
value are correspondingly high. 

^The uncertainty in the minimum and maximum hemisphere-averaged fluxes as determined from the 
two-hemisphere model flt to the phase curve (§2.4) includes the uncertainties from both the secondary 
eclipse depth estimate and the relative increase/decrease in flux from the secondary eclipse to the maxi- 
mum/minimum hemisphere-averaged flux value; the values quoted in the text of §2.4 do not include the 
additional uncertainty from the secondary eclipse depth estimate. 



series affected by the residual detector ramp. In our pre- 
vious observations of HD 189733b at 8 /im (Kinitson 
et al. 2007), long-term trends in the time series were 
caused by low-illumination pixels at the edges of our 
aperture. These pixels converged to a constant value 
more slowly than the high-illumination pixels at the cen- 
ter of the star's psf, a behavior that could be observed by 
comparing the slope of the ramp for smaller and larger 
apertures (because larger apertures contained more low- 
illumination pixels, the time series using these apertures 
will have a significant ramp even at late times). In our 
new prcflashcd observations of HD 149026 the behavior 
of the residual ramp in the data does not appear to be 
correlated with the illumination level of individual pix- 
els, as we see identical behaviors for apertures ranging 
from 3.5 — 7.0 pixels once the background flux has been 
subtracted. We take this as evidence that the low-flux, 
large-slope component of the ramp has been effectively 
removed from the data, and proceed to fit the data start- 
ing 10 hours into the observations without any further 
corrections. 

We first fit a simple linear function of time to the data, 
taking all points between 10 < i < 36.53 hours after the 
start of the observations and excluding all data beginning 
8 minutes {2a uncertainty in the original best-fit eclipse 
time, see §3.1.2) before the start of ingress and ending 8 
minutes after the end of egress in the best-fit secondary 
eclipse light curve. This corresponds to a range in or- 
bital phase of 0.085 — 0.534. We find an increase in flux 
of 0.0227% ± 0.0066% over this interval, with a signifi- 
cance of 3.4fT, for our best-fit linear solution, which has 
a value of 228,872.69, two degress of freedom, and 
228,920 points in the fit. Moving to a quadratic func- 



tion decreases the value of the fit by 4.3, indicating 
that there is detectiblc curvature in the shape of the ob- 
served phase variation. For the quadratic fit we find a 
flux minimum located at an orbital phase of 0.20 ± 0.17, 
corresponding to a central meridian longitude 72 ± 61° 
west of the antistellar point. 

We also fit a more realistic two-hemisphere model 
(Cowan & Agol 2008) to these same data. We first 
fit a model where we fix the locations of the day and 
night side hemispheres to be centered on the substel- 
lar and antistellar points on the planet, respectively. 
Figure 7 shows the binned data and resulting best-fit 
phase curve (dotted line). We find a maximum (day- 
side) hemisphere-averaged flux of 0.0305% ± 0.0031% 
and a minimum (nightside) hemisphere- averaged flux of 
0.0118% ± 0.0036% from this fit. Note that this model 
underestimates the flux at secondary eclipse (this is anal- 
ogous to the maximum hemisphere-averaged flux in this 
case) in order to provide a better fit to the shape of the 
phase curve prior to the start of the secondary eclipse. 
The value of this fit is marginally worse than either 
the linear or quadratic fits {5x^ = 1.2 as compared to the 
linear fit, 6x^ = 5.5 for the quadratic fit). We also tried 
fitting a model in which the day-night fiux distribution 
varied as the cosine of the longitude from the substel- 
lar point; this is the simplest version of the sinusoidal 
class of models described in Cowan & Agol (2008). The 
resulting light curve was almost identical to that of the 
two-hemisphere model, with 5x^ = 0.03 (a slightly worse 
fit). 

We obtain an improved fit {Sx^ = —5.1 as compared 

to the linear fit and (5x^ = —0.9 as compared to the 
quadratic fit) if we allow the locations of the two hemi- 
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spheres to vary in our fit (solid line in Fig. 7). In this 
fit, the maximum and minimum hemisphere- averaged 
fluxes do not necessarily occur at the points of supe- 
rior and inferior conjunction, although these two events 
will always occur exactly half an orbit apart for the case 
where the planet has a circular orbit. The maximum 
and minimum hemisphere-averaged fluxes in this fit are 
0.0539% ±0.0115% and 0.0161% ±0.0024%, respectively. 
In this model the centers of the two hemispheres are 
shifted to central meridian longitudes 73 ±17° west of the 
substellar and antistellar points. The artificially small 
uncertainty on this phase shift comes from the assump- 
tion of a specific physical model (namely two hemispheres 
of constant surface brightness, producing a phase curve 
where the maximum fltix always occurs 180° after the 
minimum flux). If we relax this restriction by increas- 
ing the number of slices in the fit to cither 3 (120°-wide 
slices) or 4 (90°-wide slices), we find the faintest slices 
are located at central meridian longitudes 72 ± 49° and 
88 ± 40° west of the antistellar point, respectively. This 
suggests that the true uncertainty in the location of the 
darkest region on the planet should be 40 — 50° in lon- 
gitude, comparable to and slightly smaller than the un- 
certainty from our original quadratic fit. Although we 
obtain consistent values for this location from all of our 
fits, we conclude that the offset from the antistellar point 
is not statistically significant. We also repeat the same 
fit with a model in which the flux distribution varies as 
the cosine of the longitude from the substellar point plus 
a phase offset; as before, the resulting fit was indistin- 
guishable from the best-fit phase-shifted two-hemisphere 
model {Sx"^ = -0.01). 

In our discussions below we use the linear fit to esti- 
mate the increase in flux from the planet beginning 10 
hours after the start of our observations and ending af- 
ter the end of the secondary eclipse and the correspond- 
ing ratio of the minimum to the maximum hemisphere- 
averaged fluxes. Although the two-hemisphere phase- 
shifted model provides a better fit to the data, its esti- 
mate of the maximum hemisphere-averaged flux involves 
an extrapolation to times after the end of our observa- 
tions. The linear flt makes no such assumptions, and 
allows us to make a simple statement about the increase 
in flux observed during the fraction of the planet's orbit 
spanned by our observations. 

2.4.1. Variability From Starspots 

As discussed in Knutson et al. (2007, 2009a), changes 
in the brightness of the planet-hosting star HD 189733 
as starspots rotate in and out of view can contribute to 
the measured changes in flux over the 30 — 40 hours of 
a typical phase curve observation. Although the ampli- 
tude of these variations is reduced at mid-IR wavelengths 
as compared to visible light, the signal from the planet 
is small enough that starspot variations cannot be ne- 
glected. HD 189733 is more active than most planet- 
hosting stars discovered to date, with a rotation period 
of 13 days and a corresponding 3% flux variation in the 
Stromgren {h+y)/2 band (Winn et al. 2007b). In Knut- 
son et al. (2009a) we found that the star contributed 
approximately 1/3 and l/5th of the total measured in- 
crease in flux over half of the planet's orbit in the 8 and 
24 jjLui IRAC bands, respectively. 

We characterize the behavior of HD 149206 during the 
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FlK. 8. — Relative flux variation for HD 149026 observed with 
the Til 0.8 111 APT in the Str0mgrcn y band (filled circles) around 
the time of our 8 /iiii Spitzer phase curve observations (gray shaded 
region). We assume an uncertainty of 0.0009 mag per point (0.08%, 
calculated as tfie standard deviation of the fluxes measured in a 
single night and taking the average of these values for all nights 
with four or more observations) and fit a linear function of time 
to these data. Using this fit, we derive a 2<j upper limit of 0.005% 
on the potential increase in the star's flux in Str0mgren y during 
our phase curve observations, corresponding to a 2cr upper limit of 
0.0006% at 8 jum. This is insignificant relative to the amplitude of 
the observed phase variation. 

epoch of our phase curve observations using the same 
Stromgren {h+y)/2 band as for HD 189733. The 564 ob- 
servations in our analysis were acquired between UT 2005 
April 25 and UT 2008 June 7 from an ongoing monitoring 
program carried out with the Til 0.8 m automatic photo- 
metric telescope (APT) at Fairborn Observatory (Henry 
1999; Eaton et al. 2003). Differential brightness mea- 
surements were made as the telescope nodded between 
HD 149026 and three comparison stars of comparable 
or greater brightness, following the methods described 
in Henry & Winn (2008). These data were used previ- 
ously by Nutzman et al. (2009), who computed the peri- 
odogram of the time series for periods between 0.5 — 100 
days and found no evidence for any significant periodic- 
ities. These authors placed an upper limit of 0.001 mag 
(0.1%) on the peak-to-peak amplitude of any sinusiodal 
flux variations in this period range. 

The Stromgren [b+y) /2 APT data for HD 149026 in- 
dicate that spot-induced variability in HD 149026 is a 
factor of 30 or more smaller than the variability observed 
in the same Stromgren {b+y)/2 band for HD 189733 
and therefore is unlikely to contribute signiflcantly to 
the measured 8 /itm flux increase of 0.0227% ± 0.0066% 
for this system. Because the APT data span the epoch 
of our 8 /Ltm phase curve observations, we derive a more 
specific constraint on potential flux variations over the 
10 days immediately before and after our phase curve 
observations (see Fig. 8). We use the Stromgren y data 
alone for this analysis, as they are closer to the infrared 
wavelengths of our Spitzer photometry than the averaged 
Stromgren (b+y) / 2 data. We set the relative uncertainty 
on each y band point equal to 0.0009 mag (0.08%, cal- 
culated as the standard deviation of the fluxes measured 
in a single night and taking the average of these val- 
ues for all nights with four or more observations). We 
then fit a linear function of time to the APT y band 
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data and find that the la upper hmit on the increase in 
HD 149026's flux over the 31 hours of the Spitzer obser- 
vations is 0.005% in Stromgren y. Using the conversion 
factor we derived in Knutson et al. (2008), we scale the 
upper limit y band variation to the 8 /itm band and find a 
maximum possible 8 /xm stellar flux increase of 0.0006%, 
a mere 3% of the observed rise in flux in the 8 fim band. 
As a result we ignore stellar variability for the remain- 
der of this analysis and assume that any changes in the 
measured 8 /xm flux are due to the planet and not the 
star. 

3. DISCUSSION 

3.1. Eclipse Fits and Evidence for Variability 

As discussed in the introduction, we compare our tran- 
sit and secondary eclipse fits to equivalent 8 /xm observa- 
tions of the transit of HD 149026 (Nutzman et al. 2009, 
obtained in 2007) and the secondary eclipse (H07, ob- 
tained in 2005). We use these older data to search for 
changes in either the measured radius or the 8 fj.m day- 
side flux from the planet, and to derive improved esti- 
mates for these parameters where possible. 

3.1.1. Transit and Secondary Eclipse Depth 
Comparisons 

We flnd a best-fit planet-star radius ratio of 0.05216 ± 
0.00078 when all parameters in our fit are allowed to vary 
freely (i.e., no assumptions about a/i?* or i). This is 
within 0.5(7 of the value given by Nutzman et al. (2009). 
Similarly, when wc fix a/i?,* and i to their best-fit values 
from Carter et al. (2009), we find a planet-star radius ra- 
tio of 0.05253 ± 0.00076, 0.6ct away from the 8 jim value 
given by Carter et al. (2009). We can reduce the un- 
certainties in our estimate of Rp/R^ by combining the 
values from our fit and the 8 fim transit presented in 
Nutzman et al. (2009) and re-analyzed by Carter et al. 
(2009); taking the error- weighted average of our value 
and the updated 8 /im value from Carter et al. (2009), 
we find Rp/R^ = 0.05224±0.00056. Our results are con- 
sistent with those of Carter et al. (2009), and the smaller 
error bars provided by our new measurement provides 
additional support for the trend noted by Carter et al. 
(2009), in which the transit depth between 1 — 2 /im ap- 
pears to be marginally deeper than the equivalent values 
for visible-light photometry (Winn ct al. 2008) and the 
Spitzer 8 fxm bandpass. Carter et al. (2009) conclude 
that this discrepancy is most likely due to systematic 
errors rather than a real variation in the properties of 
the planet's atmosphere; the noise in their 1.4 /im NIC- 
MOS data is a factor of two above the predicted photon 
noise and its properties are poorly understood. The good 
agreement between our 8 /zm transit and that of Nutz- 
man et al. (2009) suggest that the problem most likely 
lies with the NICMOS data, as the noise in the Spitzer 
data sets is only a factor of 1.28 and 1.15 above the pho- 
ton noise limit, respectively. 

If we compare our new secondary eclipse depth to that 
reported by H07 in this bandpass, it initially appears as 
if this quantity might be varying in time. Our best-fit 
secondary eclipse depth is 0.0411%±0.0076%, a factor of 
two smaller (3.0(j) than the value of 0.082%;° °'|'^^° from 
H07. Fixing the time of secondary eclipse in our fits to 
its predicted value results in a slightly shallower eclipse 



depth, with a net change of 0.006% or 0.8a" as compared 
to our original fit. 

It is possible that the removal of large-amplitude de- 
tector effects in the older H07 data might have resulted 
in an over-estimate of the true eclipse depth, in which 
case our two measurements might in fact be consistent. 
Although H07 used the same 0.4 s exposure time and 
8 /im IRAC subarray as wc did, our observations place 
the star at a single pixel at the center of the subarray 
whereas the H07 data cycle between nine different nod 
positions, with a total of twelve complete cycles over 6.4 
hours of observations. Each nod position has a slightly 
different gain and the move to a new nod position results 
in a slight reset of the detector ramp effect described in 
§2.2, therefore H07 must fit a total of 23 free parame- 
ters in their fit as compared to our four free parameters. 
They must also remove a ramp spanning the entire data 
set with a total rise of 1.3%, fifteen times larger than 
their inferred 0.084% decrease in flux during secondary 
eclipse. 

As compared to the H07 flts, our analysis is much sim- 
pler. We do not nod, therefore wc have no need to nor- 
malize to different nod positions. We preflash our data 
to mitigate the ramp and the secondary eclipse occurs at 
the end of our observations, where the slope is minimal 
(the best-flt coefficient for the linear function of time in 
our eclipse fits differs from zero by only 0.2(7 and we ob- 
tain the same results when we limit our fits to a simple 
constant term to describe the out-of-cclipsc flux). We 
also have more data out of eclipse than H07, which is 
helpful for determining a robust out-of-eclipse baseline. 
Although the uncertainties in our measured eclipse depth 
are only slightly smaller than those reported by H07 for 
the 2005 eclipse, wc conchide that our new eclipse depth 
is a more reliable estimate of the planet's 8 /xm flux for 
the reasons listed above. 

As discussed in §2.3.1, wc evaluate the robustness of 
the eclipse depth provided by H07 by downloading the 
images from 2005 from the Spitzer archives and perform- 
ing our own independent analysis of these data. Depend- 
ing on our specific choices for which data to trim and 
what functional forms to fit, we find eclipse depths rang- 
ing from 0.05 — 0.08%, although an eclipse depth around 
0.05% is preferred in a majority of our fits (Fig. 6). If we 
repeat these fits using the same photometry as H07 we 
find solutions with eclipse depths closer to 0.05% where 
the value is indistinguishable from the best value 
obtained by H07. Our results suggest that the 2005 data 
is in fact consistent with our new 2008 obsc;rvations, and 
we therefore conclude that there is no evidence for vari- 
ability in the planet's 8 /xm flux. 

3.1.2. Transit and Secondary Eclipse Timing 

Our best-flt transit time occurs 1.18 ± 0.62 minutes 
earlier than predicted using the ephemeris from Carter 
et al. (2009). We combine this new transit value with 

previous transit times from Winn et al. (2008), Nutz- 
man et al. (2009), and Carter et al. (2009) to derive 
an improved ephemeris for the planet (Fig. 9). We 
flnd a period of 2.8758925 ± 0.0000023 days and a mid- 
transit time of 2454597.70645 ±0.00018, consistent with 
the ephemeris from Carter et al. (2009). The NICMOS 
transits described in that work have a higher timing pre- 
cision than the Spitzer transits and the authors note that 
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Fig. 9. — Measured transit times for HD 149026b along with 
our new best-fit ephemeris. The x axis gives the number of or- 
bits between observations, where our new value is used as the zero 
point. The y axis shows the difference between the observed and 
calculated transit times. Early points with larger error bars are 
ground-based data from Sato et al. (2005), Charbonncau ct al. 
(2006), and Winn et al. (2008), while the point at -94 orbits is 
the Spitzer transit from Nutzman et al. (2009), and the four points 
between —49 and —18 orbits are HST observations from Carter et 
al. (2009). Dashed lines show the Icr uncertainties in the predicted 
transit times. The three most recent transits differ from their pre- 
dicted values by — 3.0<t, 2. Oct, and 1.6ct, respectively, if we neglect 
the correlated uncertainties in the predicted transit times. 

their measured transit times are marginally inconsistent 
(X^ = 20.16 with 14 degrees of freedom) with a con- 
stant period. In particular, the two most recent NIC- 
MOS transits, which occur 92 and 52 days prior to our 
observations, differ from their predicted values by — 2.6(t 
and -f 2.3ct (1.2 and 1.6 minutes), respectively. With our 
new ephemeris we find a value of 23.7 with 15 degrees 
of freedom, with the two most recent NICMOS transits 
and our new transit differing from their predicted values 
by — 3.0(7, 2.0(7, and 1.6(7, respectively, where we have 
neglected the correlated uncertainty in the ephemeris. 

If HD 149026b has a circular orbit, the secondary 
eclipse should occur half an orbit after the transit, with 
an additional 45 s delay from the increased light-travel 
time in the system (Loeb 2005). HOT find that their 
8 /um secondary eclipse occurred — 3.0t2'.5 minutes from 
the predicted time (i.e., earlier than predicted). Our new 
8 fxm secondary eclipse, which was observed three years 
after that of H07, differs from the predicted time in our 
nominal fits by —20.9'^\'q minutes with a significance of 
4.9(7. Because transit and secondary eclipse times are 
particularly sensitive to the presence of correlated noise 
in the data, we use the "prayer bead" method (Moutou 
et al. 2004) to calculate a second, independent estimate 
for the uncertainty in our best-fit secondary eclipse time. 
First we take the residuals from our best-fit secondary 
eclipse solution and shift them forward by ten points at 
a time, with the data at the end of the time scries wrap- 
ping around to the beginning. After each shift we add 
the best-fit eclipse function bac;k in and then fit the re- 
sulting light curve to determine the best-fit secondary 
eclipse depth and time. We repeat this process until the 
residuals have cycled around to their original positions, 
and then plot a histogram of the resulting eclipse times. 
The resulting distribution is close to Gaussian, consis- 
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Fig. 10. — Probability distribution of timing offsets from a series 
of MCMC fits to the secondary eclipse plotted in Fig. 5. These 
fits exclude the first 10 (top left), llth-20th (top right), 21st-30th 
(middle left), 31st-40th (middle right), 41st-50th (bottom left), 
and 51st-60th frames (bottom right) of each set of 64 images. The 
differing distributions of timing offsets reflect the presence of corre- 
lated noise in the data; the best-fit eclipse depth and corresponding 
uncertainty, which are less sensitive to correlated noise on shorter 
time scales, are the same in all of these fits. We adopt the largest 
timing errors from this set of fits (middle left histogram, -1-7.1 and 
—6.4 minutes) as a conservative estimate of the uncertainty in the 
best-fit time of secondary eclipse (see discussion in §3.1.2). For our 
best-fit eclipse time we keep the original value from our fit to the 
full data set. 

tent with the results of our MCMC analysis, and has a 
standard deviation of 5.0 minutes. 

As discussed in §2.1, subarray photometry is recorded 
in sets of 64 images, and the first 5 — 10 frames and 58th 
frames in each set of 64 tend to have low flux values. 
In our initial analysis we concluded that this effect was 
removed by the background subtraction, but it is pos- 
sible that some correlated noise remains on these time 
scales. As a test we re-ran otir secondary eclipse MCMC 
fits without the first 10 and 58tli frames, and found that 
we obtained identical values for the eclipse depth, but 
the distribution of values for the best-fit eclipse time was 
both broader and noticeably asymmetric, with Icr uncer- 
tainties on the best-fit eclipse time of Zl'y- In order to 
determine if this change was connected to the presence 
of the low-background frames, we re-ran the same fits 
excluding the first 10, llth-20th, 21st-30th, 31st-40th, 
41st-50th, and 51st-60th frames of each set of 64 images 
(sec Fig. 10). In all cases we obtained consistent values 
for the best-fit eclipse depth and corresponding uncer- 
tainties, as well as the best-fit eclipse time, but the un- 
certainties on this eclipse time ranged from (exclud- 
ing the 51st-60th frames) to tgj (excluding the 21st-30th 
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frames). As a more conservative error estimate, we use 
the largest of these error bars for the best-fit ecHpse time 
in our subsequent analysis below. This uncertainty sug- 
gests that the measured eclipse time is only marginally 
inconsistent (2.9ct as opposed to our earlier value of 4.9cr) 
with the predicted value for a circular orbit. The uncer- 
tainty in the secondary eclipse depth is unchanged, as 
we obtain consistent results in all cases. We also re- 
peated this same test for the transit fits and found that 
the uncertainties in the best-fit planet-star radius ratio 
and transit time were similarly consistent. We conclude 
that the only parameter sensitive to this effect is the sec- 
ondary eclipse time, a fact that is reflected in our new, 
larger error estimate for this quantity. 

Using these new error bars, we find that the offset in 
our best-fit secondary eclipse time differs from that of 
H07 by 2.3ct. In our re-analysis of the H07 data (§3.1.1) 
we find a best-fit eclipse time 3.3 minutes later than the 
value derived by H07, albeit with larger uncertainties 
(±3.5 minutes as opposed to H07's value of j^2'g minutes). 
Using this value, the timing offset in the 2005 eclipse 
differs from that of our 2008 eclipse by 2.6a. If this offset 
is confirmed by future observations, it would suggest that 
the planet's orbit may have changed during the three year 
interval between the two observations. 

The timing offset observed in our new data corresponds 
to a value of — 0.0079to;oo25 for e cos (w) where e is the 
planet's orbital eccentricity and lo is the argument of 
pericenter. This is consistent with the radial velocity 
observations of this system from Sato et al. (2005) and 
Wolf et al. (2007), which are poorly sampled and thus 
do not provide strong constraints on the planet's orbital 
eccentricity. Madhusudhan & Winn (2009) derive upper 
limits on e sin (w) and e cos (w) from a simultaneous fit to 
the available transit, secondary eclipse, and radial veloc- 
ity data, but the uncertainty in e cos (to) is dominated by 
the uncertainty in the secondary eclipse time from H07 
and their value for esin(tj) = 0.109to;o68 is still poorly 
constrained relative to our estimate for e cos (u) . If real, 
the observed change in e cos (cu) could be due to either 
a changing orbital eccentricity or a precessing orbit (i.e., 
a change in e or in w). The addition of four new sec- 
ondary eclipses observed as part of an ongoing Spitzer 
target of opportunity program (J. Harrington 2009, pri- 
vate communication) should either confirm or disprove 
the presence of the timing offset suggest by our data, 
as well as constraining any temporal variability in this 
quantity. 

3.2. Comparisons to ID Atmosphere Models 

Spitzer observations of secondary eclipses for close-in 
hot Jupiters (e.g. Deming et al. 2005, 2006; Charbon- 
neau et al. 2005, 2008; Harrington et al. 2007; Richard- 
son et al. 2007; Grillmair et al. 2007, 2009; Machalek et 
al. 2008; Knutson et al. 2008, 2009b) suggest that their 
atmospheres can be divided into two classes, depending 
on whether or not they have a temperature inversion and 
water emission (as opposed to water absorption) features 
in the near- to mid-IR. Planets with inversions appear 
to require an additional absorber at altitude, preferably 
with a large cross-section in visible light where the stars 
emit most of their flux, in order to create the inversion. 
Several authors (Hubcny et al. 2003; Fortney et al. 2006, 
2008; Burrows et al. 2007b, 2008a) have suggested that 
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Fig. 11. — Best-fit 8 fim secondary eclipse depths from this paper 
(filled circle) and from Harrington et al. (2007) (filled star) with 
ztlcr error bars overplotted. We re-analyze the Harrington et al. 
(2007) data in this paper and find a smaller value consistent with 
our new secondary eclipse depth. Two dayside atmosphere mod- 
els for HD 149026b are shown for comparison (these are similar 
to the models in Fortney et al. 2006, but these models have been 
updated to match the revised system parameters from Carter et 
al. 2009). The black model (solid line) assumes half the incident 
flux is transported to the night side (i.e. no day-night temperature 
gradient), does not include TiO opacity and therefore has no in- 
version, and shows water and CO bands in absorption. The other 
black model (dotted line) makes the same assumptions but uses 
a 30x solar metallicity atmosphere for the planet instead of solar 
metallicity. This assumption has a minimal effect on the resulting 
planet-star flux ratio. In contrast to the two black models, the 
red model absorbs and re-radiates all of the incident fiux from its 
star on the day side (this assumes effectively zero fiux from the 
cold night side), includes TiO opacity, and has formed a tempera- 
ture inversion with water and CO features in emission. The filled 
squares give the predicted planet-star fiux ratios integrated over 
the six Spitzer bandpasses (3.6, 4.5, 5.8, 8.0, 16, and 24 ^m) for 
each of these models. 



gas-phase TiO might provide the required opacity, al- 
though cold trap effects would reduce the amount of TiO 
in the upper atmosphere to much lower levels unless there 
is fairly vigorous mixing in this region (Spiegel et al. 2009; 
Showman et al. 2009). In this picture planets with suffi- 
ciently high levels of incident flux would have gas-phase 
TiO and would form inversions, whereas cooler planets 
would not. Although this theory works well for two of the 
three hot Jupiters with confirmed inversions (Knutson et 
al. 2008, 2009b), it breaks down for XO-lb (Machalek et 
al. 2008), which is too cool for gas-phase TiO and yet ap- 
pears to have an inversion. Zahnle et al. (2009) recently 
proposed an alternative theory where photochemistry in 
the upper atmosphere might lead to the production of 
enough disequilibrium sulfur compounds to form an in- 
version. In this case the authors find that the formation 
rates for the sulfur compounds of interest appear to be 
independent of stellar insolation, atmospheric metallic- 
ity, and temperature in the 1200 — 2000 K range spanned 
by their models. Although these models correctly predict 
that XO-lb should have an inversion, they fail to explain 
why HD 189733b and TrES-1 do not. 

HD 149026b is a clear candidate for an inversion in ei- 
ther scenario, as it receives enough incident flux in the 
zero-albedo case to have gas-phase TiO in the hotter re- 
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gions of its dayside atmosphere, and the sulfur photo- 
chemistry model suggests that all close-in planets should 
have inversions. Although this planet has an equilib- 
rium temperature of 1740 K, the 8 fxm secondary eclipse 
depth measured by H07 corresponds to an 8 /xm bright- 
ness temperature of 2300 ± 200 K. In order to match this 
high dayside flux, models must keep as much energy as 
possible on the day side (i.e., a large day-night tempera- 
ture gradient). Models with a temperature inversion and 
water in emission also provide a better match to the data 
(Fortney et al. 2006; Burrows et al. 2008a), as they tend 
to have higher 8 /im fluxes than their non-inverted coun- 
terparts, which have water absorption features in this 
wavelength range (see Fig. 11). With only one wave- 
length of data, however, we cannot directly verify the 
presence or absence of a temperature inversion in HD 
149026b's dayside atmosphere. 

We measure a considerably smaller 8 /xm secondary 
eclipse depth in our new 8 /itm observations; our value as 
plotted in Fig. 11 is half that reported by HOT, with a 
brightness temperature of only 1440 ± 150 K. This value 
is best matched by models with strong circulation be- 
tween the two hemispheres (i.e., half the incident flux 
is transported to the night side), as this serves to re- 
duce the predicted dayside temperatures. Models with 
absorption features in the 8 /xm bandpass, such as the 
non- inverted model in Fig. 11, also tend to provide a 
better fit to the data, but the observed 8 /xm flux could 
also be matched by a model with blackbody emission and 
a non-zero albedo. It is risky to extrapolate too much 
on the basis of observations at a single wavelength, and 
definite confirmation of the presence or absence of an in- 
version must await the addition of observations at other 
wavelengths. Wc will carry out full-orbit phase curve 
observations of this system in the 3.6 and 4.5 /xm IRAC 
bands, which have proven to be particularly useful for 
diagnosing the presence of such temperature inversions 
(e.g. Knutson et al. 2008), as part of Spitzer's extended 
warm mission from 2009 — 2011. 

Although our reanalysis of the HOT data is consis- 
tent with a constant 8 /tm secondary eclipse depth from 
2005 — 2008, it is interesting to speculate on the possi- 
bility that HD 149026b is time variable. Very few plan- 
ets have more than one secondary eclipse observation in 
the same IRAC channel, but this planet has been ob- 
served four times at 8 /tm and twice at 5.8 /tm (two of 
the 8 /xm secondary eclipses and two new 5.8 /xm sec- 
ondary eclipses were observed as part of an ongoing 
Spitzer target of opportunity program and are currently 
unpublished; J. Harrington 2009, private communica- 
tion). These new data, when combined with our eclipse 
and that of HOT, will allow for a search for variability on 
a time scale of four years with multiple epochs of obser- 
vations, and will be sensitive to smaller variations than 
those detectible in the current data. 

Time variability is common in planetary atmospheres, 
and several hot Jupiter models have suggested highly 
time-variable flows (e.g., Cho et al. 2003, 2008; Lang- 
ton & Laughhn 2008; Menou & Raushcr 2009). Obtain- 
ing large-amplitude time variability on the global scale 
requires the existence of dynamical instabilities that (i) 
have length scales comparable to (or larger than) a plan- 
etary radius, and (ii) having instability growth times 
shorter than the timescales of any processes that might 



damp the instabilities; in practice this probably requires 
growth times shorter than the relevant radiative time 
constants. To date, most theoretical studies suggesting 
that global-scale variability could be significant adopt 
two-dimensional models; some of these models essentially 
ensure a time-variable flow by adopting a turbulent ini- 
tial condition (Cho et al. 2003, 2008; Rauschcr ct al. 
2008) while others generate highly fluctuating flows even 
when initialized from rest (Langton & Laughlin 2008). 
On the other hand, most three-dimensional models of 
synchronously rotating hot Jupiters on circular orbits 
published to date produce relatively steady flow patterns 
that do not exhibit large-amplitude global-scale variabil- 
ity (Showman & Guillot 2002; Cooper & Showman 2005, 
2006; Showman et al. 2008, 2009; Dobbs-Dixon & Lin 
2008). One three-dimensional study produced highly 
time variable flows (Menou & Raushcr 2009); the differ- 
ent behavior might result from the very different forcing 
scheme adopted in Menou & Rauscher (and their place- 
ment of the bottom surface at relatively shallow pressure) 
as compared to previous 3D studies. These theoretical 
predictions for variability will likely evolve as models ex- 
plore a wider parameter space in forcing schemes and 
planetary parameters, while the availability of multiple 
epochs of secondary eclipse observations for a few select 
planets should provide the first observational tests for 
such variability. In Agol et al. (2009) we present one 
such search for variability in the 8 /xm emission from 
HD 189733b, where we conclude that this planet's emis- 
sion appears to be constant at the level of 10% over the 
five secondary eclipses observed. 

3.3. The 3D Picture: A Global Energy Budget 

Rather than extrapolating the nightside temperature 
based on the ID dayside models plotted in Figure 11, we 
can use our measurement of the planet's phase variation 
to directly characterize this quantity. We measure an in- 
crease in flux of 0.022T% ± 0.0066% over the course of 
our observations; when combined with our new value for 
the secondary eclipse depth, which gives the dayside flux, 
this corresponds to a minimum hemisphere-averaged flux 
that is 45% ± 19% of the maximum hemisphere-averaged 
flux. If we assume the planet emits as a blackbody and 
solve for the temperature needed to match the maximum 
and minimum hemisphere-averaged planet-star flux ra- 
tios (Table 2), this corresponds to a maximum (dayside) 
hemisphere-averaged brightness temperature of 1440 ± 
150 K and a minimum (nightside) hemisphere- averaged 
brightness temperature of 960 ± 240 K. The day- night 
temperature gradient in this case is then 480 ± 140 K, 
where the uncertainty in this quantity simply reflects the 
uncertainty in the measured rise in flux during these 
observations. The uncertainty in the dayside temper- 
ature is derived from the uncertainty in the depth of 
the secondary eclipse, while the uncertainty in the night- 
side temperature is a function of both the uncertainty 
in the secondary eclipse depth and the day-night flux 
difference and is therefore larger than cither of these 
quantities. The resulting day-night temperature gradi- 
ent for HD 149026b is larger than the approximately 
200 ± 50 K gradient observed for HD 189T33b (Knut- 
son et al. 2009a), but this difi'erence is not statistically 
significant. 

If we instead use the two-hemisphere model fit to the 
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Fig. 12. — This plot shows the best-fit 8 fim night-side flux 
from our phase curve flt with ztlcr error bars overplotted. The two 
overplotted models show the predicted emission from an object in 
isolation (i.e., zero incident flux) with an effective temperature of 
1300 K (black line) or 1600 K (red line). The dotted black line is 
the best-flt dayside model from Fig. 11. The flUed squares give 
the predicted planet-star flux ratios integrated over the six Spitzer 
bandpasses (3.6, 4.5, 5.8, 8.0, 16, and 24 fira) for each of these 
models. 

phase curve from §2.4 to determine the maximum and 
minimum hemisphere-averaged fluxes, we find a larger 
gradient of 780 ± 240 K between these two hemispheres. 
The maximum flux in this model occurs after the end 
of our observations (see Fig. 7), and the increased tem- 
perature gradient is a direct result of this extrapolated 
maximum flux value. We prefer the linear fit as a more 
conservative estimate of the day-night temperature gra- 
dient, although it is in some sense a lower limit as the 
planet still appeared to be increasing in brightness at the 
end of our observations. 

Although both the minimum and maximum 
hemisphere-averaged temperatures are lower than 
HD 149026b's 1770 K equilibrium temperature assum- 
ing an albedo of zero and uniform (47r) re-radiation of 
incident flux, this does not necessarily mean that there 
is a problem with its energy budget. If the planet has an 
emission spectrum with water absorption it will appear 
correspondingly fainter in our 8 /im bandpass, and in 
fact we see that the ID model spectrum in Fig. 11 is a 
good match for our dayside planet-star flux ratio. This 
model has a Bond albedo of near zero and an effective 
temperature of 1750 K, and it assumes that half the 
incident flux is redistributed to the night side. We 
obtain comparable fits with both solar-metallicity and 
30x solar metallicity models for the planet's atmosphere. 

For the planet's night side, which receives no incident 
flux from the star, we use a model for an object in isola- 
tion and vary the effective temperature in order to match 
the observed 8 /im flux (see Fig. 12). We find that a 
model with Teff — 1300 K provides the best fit, although 
Teff — 1600 K is still within the la uncertainty for the 
night-side fiux. Using an effective temperature of 1750 K 
for the day side and 1300 K for the night side we find that 
the planet emits 60% as much fiux as it receives from its 
parent star in the zero-albedo case. In order to balance 
these two quantities the planet must then have a Bond 



albedo of 40%. If we carry out the same calculation with 
a 1600 K night side we find a Bond albedo of 17%. This 
is not a fully self-consistent solution, as we do not re-run 
our dayside model with a non-zero albedo (which could 
perhaps come from refiecting clouds) and the appropri- 
ate redistribution required to match the nightside flux, 
but the difference is insignificant in light of the current 
uncertainties in the measured day- and night-side fluxes 
and the limited wavelength range of the available data. 

Upper limits on reflected light from hot Jupiters such 
as HD 209458b (Rowe et al. 2008) suggest that their albe- 
dos are small (the la upper limit for the geometric albedo 
of HD 209458b is 0.08), a conclusion that is in good 
agreement with the predictions of atmosphere models for 
these planets (e.g. Burrows et al. 2008b). CoRoT-lb was 
also recently found to have a geometric albedo < 0.20 
(this planet still has some thermal emission in the red 
CoRoT band of these observations, and the observed sec- 
ondary eclipse signal could be due to either reflected light 
or thermal emission from the planet; Snellen et al. 2009), 
placing it in a regime similar to that of HD 209458b. 
HD 149026b appears to be an exception to this trend, 
although it is worth remembering that the planet emits 
most of its flux at shorter wavelengths where it may be 
brighter than predicted by our atmosphere models. If 
this was the case, the planet's emission could still be 
consistent with an albedo of zero. Additional observa- 
tions at shorter wavelengths should allow us to test the 
validity of our atmosphere models and to provide tighter 
constraints on the planet's Bond albedo. 

4. CONCLUSIONS 

In this work we present a revised value for the 
8 /im dayside flux from the hot Saturn HD 149026b and 
the flrst estimate of the nightside flux and correspond- 
ing day- night temperature gradient. Our new secondary 
eclipse depth is half that reported by H07 in the same 
8 /im bandpass for observations obtained in 2005, and 
is best matched by atmosphere models with strong day- 
night circulation, water absorption in the 8 //m band- 
pass, a non-zero albedo, or some combination of the 
three. Our conclusion stands in direct contrast to that 
reached by H07, who concluded that the planet's appar- 
ent high flux in the 8 fim bandpass was best described 
by a model with a large day-night temperature gradi- 
ent and/or water emission features. Our re-analysis of 
the H07 data indicates that the 2005 data are in fact 
consistent with our new, smaller 8 /im eclipse depth, 
and we conclude that there is no evidence for variability 
in the planet's 8 /xm flux over the three years between 
these two observations. There have been two additional 
8 /im secondary eclipse observations (currently unpub- 
lished) obtained within the past year; these data should 
provide a more rigorous test for any potential variability 
in HD 149026b's dayside emission. Observations of the 
secondary eclipse at 3.6 and 4.5 /im, which have been 
scheduled as part of Spitzer 's extended two-year warm 
mission, will provide an important complement to these 
data and should allow us to conflrm the presence or ab- 
sence of a temperature inversion. 

Although we find no evidence for changes in the 
planet's dayside flux, our new secondary eclipse oc- 
curs — 20.9lg'5 minutes earlier than the predicted time, 
2.3a earlier than the 2005 eclipse from H07 (the best- 
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fit eclipse time from our re-analysis of the 2005 data is 
comparable to that reported by H07). The timing off- 
set for our new secondary eclipse corresponds to a value 
of — 0.0079lo!oo25 fo'' ecos(w) where e is the planet's or- 
bital eccentricity and u> is the argument of pericenter. 
By mid-2011 Spitzer will have observed eight secondary 
eclipses for this planet (one each at 3.6 and 4.5 /im, two 
at 5.8 jim, and four at 8.0 y^m); these data will allow us 
to confirm or disprove the presence of an offset in the 
secondary eclipse times and to search for any evidence of 
temporal variability. 

Our measurement of the planet's phase curve indicates 
that the minimum hemisphere-averaged 8 fim flux for 
this planet is 45% ± 19% of the maximum hemisphere- 
averaged flux. We use our estimates of the day- and 
night-side fluxes to construct an energy budget for the 
planet and conclude that it appears to emit less flux 
than it absorbs from its parent star, corresponding to 
a non-zero Bond albedo, although the models are poorly 
constrained by the current data. Because we only have 
measurements of the day- and night-side fluxes in a sin- 
gle bandpass located on the long-wavelength tail of the 
planet's emission spectrum, our observations may still 
be consistent with a zero albedo if the planet is brighter 
than predicted at shorter wavelengths. Although our 
models have similar predictions for solar and 30 x solar 
atmospheric metallicities, variations in the abundances 
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of specific elements might still affect our predictions for 
the 8 /im bandpass. 

Over the next two years we will obtain full-orbit phase 
curve observations of HD 149026b at 3.6 and 4.5 /xm as 
part of Spitzer 's extended warm mission; these wave- 
lengths are close to the peak of the planet's emission 
spectrum and will therefore provide much tighter con- 
straints on its energy budget. The full phase coverage 
will also allow us to better resolve the locations of hot 
and cold regions in the planet's atmosphere. These ob- 
servations will take place as part of a larger program to 
provide multi-wavelength characterizations of the phase 
curves of six extrasolar planets, including HD 189733b, 
HD 209458b, HD 149026b, HAT-P-2, HAT-P-7, and GJ 
436b, and comparisons between these systems will allow 
us to investigate the importance of irradiation, rotation 
rate, stirface gravity, eccentricity, and stellar metallicity 
in determining the pressure-temperature profiles and dy- 
namic meteorology of these unusual atmospheres. 
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